home *** CD-ROM | disk | FTP | other *** search
/ Aminet 1 (Walnut Creek) / Aminet - June 1993 [Walnut Creek].iso / usenet / sources / volume2 / aplictns / newton.1 < prev   
Text File  |  1988-10-29  |  20KB  |  815 lines
  1. Path: xanth!nic.MR.NET!tank!ncar!mailrus!ulowell!page
  2. From: page@swan.ulowell.edu (Bob Page)
  3. Newsgroups: comp.sources.amiga
  4. Subject: v02i037:  newton - estimate the roots of a polynomial equation
  5. Message-ID: <9924@swan.ulowell.edu>
  6. Date: 29 Oct 88 01:26:06 GMT
  7. Organization: University of Lowell, Computer Science Dept.
  8. Lines: 804
  9. Approved: page@swan.ulowell.edu
  10.  
  11. Submitted-by: barrett@cs.jhu.edu (Dan Barrett)
  12. Posting-number: Volume 2, Issue 37
  13. Archive-name: applications/newton.1
  14.  
  15. Given a polynomial equation F(x) = 0, Newton estimates the roots of
  16. the equation.  The program uses an algorithm known as "Newton's
  17. Method".  You can read about the specifics in any college textbook on
  18. Numerical Analysis.
  19.  
  20. # This is a shell archive.  Remove anything before this line
  21. # then unpack it by saving it in a file and typing "sh file"
  22. # (Files unpacked will be owned by you and have default permissions).
  23. # This archive contains the following files:
  24. #    ./Makefile
  25. #    ./Makefile.unix
  26. #    ./README
  27. #    ./complex.c
  28. #    ./decl.h
  29. #    ./ds.c
  30. #    ./io.c
  31. #    ./main.c
  32. #    ./strtok.c
  33. #
  34. if `test ! -s Makefile`
  35. then
  36. echo "writing Makefile"
  37. cat > Makefile << '\Rogue\Monster\'
  38. ############################################################################
  39. # Makefile for the "Newton" program, using Manx C V3.6a.
  40. #
  41. # Written by Daniel Barrett.  100% PUBLIC DOMAIN.
  42. ############################################################################
  43.  
  44. LFLAGS=-g +Q
  45. SRC=main.c complex.c io.c ds.c strtok.c
  46. OBJ=main.o complex.o io.o ds.o strtok.o
  47. PROG=Newton
  48.  
  49. # Use ma32.lib for IEEE double precision floating point arithmetic.
  50. CFLAGS=+L -n +Iinclude.pre +fi
  51. LIBS=-lma32 -lc32
  52.  
  53. # Or, Use m32.lib for faster, less precise arithmetic.
  54. #CFLAGS=+L -n +Iinclude.pre
  55. #LIBS=-lm32 -lc32
  56.  
  57. # Or, use 68881 library for double precision arithmetic.
  58. # Program crashed when I tried to run it after this compilation.
  59. #CFLAGS=+L -n +Iinclude.pre +f8
  60. #LIBS=-lm832 -lc32
  61.     
  62. amiga:        include.pre $(PROG)
  63.  
  64. include.pre:    decl.h
  65.         cc +L -n +Hinclude.pre decl.h
  66.         delete decl.o
  67.  
  68. $(PROG):    $(OBJ) 
  69.         ln $(LFLAGS) $(OBJ) -o $(PROG) $(LIBS)
  70.     
  71. $(OBJ):        decl.h
  72.  
  73. clean:
  74.         delete \#?.o $(PROG) $(PROG).dbg include.pre
  75. \Rogue\Monster\
  76. else
  77.   echo "will not over write Makefile"
  78. fi
  79. if [ `wc -c Makefile | awk '{printf $1}'` -ne 994 ]
  80. then
  81. echo `wc -c Makefile | awk '{print "Got " $1 ", Expected " 994}'`
  82. fi
  83. if `test ! -s Makefile.unix`
  84. then
  85. echo "writing Makefile.unix"
  86. cat > Makefile.unix << '\Rogue\Monster\'
  87. ############################################################################
  88. # Makefile for the "Newton" program, UNIX 4.2BSD.
  89. #
  90. # Written by Daniel Barrett.  100% PUBLIC DOMAIN.
  91. ############################################################################
  92. SRC=main.c complex.c io.c ds.c 
  93. OBJ=main.o complex.o io.o ds.o 
  94. LIBS=-lm
  95. PROG=newton
  96. CFLAGS=-DUNIX -g
  97.  
  98.  
  99. all:        $(PROG)
  100.  
  101.  
  102. $(PROG):    $(OBJ)
  103.         cc $(OBJ) -o $(PROG) $(LIBS)
  104.     
  105. $(OBJ):        decl.h
  106.  
  107. clean:
  108.         rm -f *.o core
  109. \Rogue\Monster\
  110. else
  111.   echo "will not over write Makefile.unix"
  112. fi
  113. if [ `wc -c Makefile.unix | awk '{printf $1}'` -ne 466 ]
  114. then
  115. echo `wc -c Makefile.unix | awk '{print "Got " $1 ", Expected " 466}'`
  116. fi
  117. if `test ! -s README`
  118. then
  119. echo "writing README"
  120. cat > README << '\Rogue\Monster\'
  121. ########################################################################
  122. # Newton:    Estimate the roots of a polynomial equation.
  123. #        Author:  Dan Barrett, barrett@cs.jhu.edu (ARPAnet).
  124. #        100% PUBLIC DOMAIN, both source code & executable.
  125. #
  126. #        Written for the Commodore Amiga, September 1988.
  127. #        Runs from the CLI only.
  128. #        Also compiles & runs under Berkeley UNIX 4.2.
  129. ########################################################################
  130.  
  131. What is Newton?
  132. ---------------
  133.     Given a polynomial equation F(x) = 0, Newton estimates the 
  134. roots of the equation.  The program uses an algorithm known as
  135. "Newton's Method".  You can read about the specifics in any college
  136. textbook on Numerical Analysis.
  137.  
  138.  
  139. How is it used?
  140. ---------------
  141.     From the CLI, type "newton", followed by <RETURN>.
  142.  
  143.     You will be prompted for the degree of the polynomial equation.
  144. The degree is the largest exponent in the equation.  The maximum
  145. allowable degree is 20.  You can change this in the source code
  146. by changing the value of MAX_DEGREE in decl.h.
  147.  
  148.     Next, you are prompted for the "desired accuracy".  Since
  149. Newton's Method is a "numerical" method, the answers you get are
  150. only estimates of the actual value.
  151.     The accuracy is a measure of how close an estimate must
  152. be before it is considered "correct".  In mathematical language,
  153. the Euclidean distance (in the complex plane) between the current
  154. estimate and the previous estimate must not exceed the "accuracy".
  155.  
  156.     In the source code, the variable "epsilon" represents the
  157. accuracy.  See also the function ErrorTooBig(), which actually tests
  158. the accuracy of the latest estimate.
  159.  
  160.     To choose the default accuracy, simply press <RETURN>.  Otherwise,
  161. type in your desired accuracy.  The smaller the value, the more accurate
  162. your result (and the longer the program will run).
  163.  
  164.     Finally, you are prompted for your polynomial coefficients.
  165. For example, when you are asked:
  166.  
  167.         X^5 coefficient?:
  168.  
  169. you should type in the coefficient of your X-to-the-5th-power term.
  170.  
  171.     Coefficients may be real or complex numbers.  If the coefficient
  172. is real, simply type it in and press <RETURN>.  If the coefficient is
  173. complex, type in its real and imaginary parts, separated by spaces
  174. and/or tabs, and press <RETURN>.
  175.  
  176. "Real" example:
  177.  
  178.         X^5 coefficient?: 6.33 <RETURN>
  179.  
  180. "Complex" example, for the complex value -5 + 0.43i :
  181.  
  182.         X^5 coefficient?: -5  0.43 <RETURN>
  183.  
  184. Note that you must type in EVERY coefficient, even if it is a zero.
  185.  
  186.     While you are entering coefficients, you can type "H", "h", or "?"
  187. for a brief "help" message.
  188.  
  189.  
  190. Now What Happens?
  191. -----------------
  192.  
  193.     First, "Newton" will display your polynomial, as you typed
  194. it in.
  195.  
  196.     Then, "Newton" will calculate the roots of your polynomial
  197. equation.
  198.  
  199.     This calculation might take some time (a few minutes), depending
  200. on the degree of the equation and the desired accuracy.  Higher accuracy
  201. takes more time, as you might guess.  Sometimes, the calculation takes
  202. only a few seconds.
  203.  
  204.     Since "Newton" uses a numerical algorithm, there is the possibility
  205. that it will not find a solution.  If "Newton" cannot find the value of
  206. a root after 10,000 iterations, it will report failure.  (This number
  207. of iterations is set in decl.h as MAX_ITERATIONS... feel free to change
  208. it.)
  209.  
  210.     So, if "Newton" seems to be sitting there not doing anything,
  211. don't panic.  After a few minutes, it will either find a solution
  212. or quit automatically.
  213.  
  214.  
  215. Compiling Information
  216. ---------------------
  217.  
  218.     The enclosed "Makefile" is for use with Aztec Manx C, V3.6a.
  219. "Newton" probably compiles & runs under Lattice C, though I can't
  220. test it.
  221.  
  222.     Math (floating point) calculations can be handled in three
  223. different ways.  Edit the Makefile and uncomment the desired
  224. CFLAGS and LIBS parameters.  They are explained in the Makefile itself.
  225. See also the "New Options For CC" page (cc.ap.1, V3.4a) in the Manx
  226. C manual.
  227.  
  228.     The enclosed "Makefile.unix" is for compiling "Newton" under
  229. UNIX (Berkeley UNIX 4.2 or Ultrix).  It does not use strtok.c; I am
  230. assuming that your UNIX has the strtok() function provided.
  231. \Rogue\Monster\
  232. else
  233.   echo "will not over write README"
  234. fi
  235. if [ `wc -c README | awk '{printf $1}'` -ne 4031 ]
  236. then
  237. echo `wc -c README | awk '{print "Got " $1 ", Expected " 4031}'`
  238. fi
  239. if `test ! -s complex.c`
  240. then
  241. echo "writing complex.c"
  242. cat > complex.c << '\Rogue\Monster\'
  243. /* complex.c:    Complex arithmetic routines.
  244.  *
  245.  * Written by Daniel Barrett.  100% PUBLIC DOMAIN. */
  246.  
  247. #include "decl.h"
  248.  
  249. /* Complex arithmetic functions Add(), Subtract(), Multiply() and Divide()
  250.  * perform as their names indicate.  They perform their operation on their
  251.  * first 2 arguments, and return the result in the third argument. */
  252.  
  253. Add(a, b, sum)
  254. complex a, b, *sum;
  255. {
  256.     sum->n[REAL] = a.n[REAL] + b.n[REAL];
  257.     sum->n[IMAG] = a.n[IMAG] + b.n[IMAG];
  258. }
  259.  
  260.     
  261. Subtract(a, b, diff)
  262. complex a, b, *diff;
  263. {
  264.     diff->n[REAL] = a.n[REAL] - b.n[REAL];
  265.     diff->n[IMAG] = a.n[IMAG] - b.n[IMAG];
  266. }
  267.  
  268.     
  269. Multiply(a, b, prod)
  270. complex a, b, *prod;
  271. {
  272.     prod->n[REAL] = (a.n[REAL] * b.n[REAL]) - (a.n[IMAG] * b.n[IMAG]);
  273.     prod->n[IMAG] = (a.n[REAL] * b.n[IMAG]) + (a.n[IMAG] * b.n[REAL]);
  274. }
  275.  
  276.     
  277. Divide(a, b, quot)
  278. complex a, b, *quot;
  279. {
  280.     double denom;
  281.  
  282.     denom = SQR(b.n[REAL]) + SQR(b.n[IMAG]);
  283.     if (denom == 0.0)
  284.         printf("Divide by zero error!\n"), exit(20);
  285.  
  286.     quot->n[REAL] = ((a.n[REAL] * b.n[REAL]) + (a.n[IMAG] * b.n[IMAG]))
  287.             / denom;
  288.     quot->n[IMAG] = ((a.n[IMAG] * b.n[REAL]) - (a.n[REAL] * b.n[IMAG]))
  289.             / denom;
  290. }
  291.  
  292.     
  293. SynthDivide(poly, comp, stop)
  294. /* Computes the synthetic division of the polynomial poly[] by (z-comp), 
  295.  * where "z" is assumed the complex variable in our polynomial. */
  296. complex poly[], comp;
  297. int stop;
  298. {
  299.     int i;
  300.     complex prod;
  301.  
  302.     for (i=1; i<=stop; i++) {
  303.         Multiply(poly[i-1], comp, &prod);
  304.         Add(poly[i], prod, &poly[i]);
  305.     }
  306. }
  307.  
  308.     
  309. Iterate(poly, zOld, n, zNew)
  310. /* One iteration of Newton's method.  "zOld" is the old guess of the value
  311.  * of a root of poly[].  "zNew" is the newly calculated guess. */
  312. complex poly[], zOld;
  313. int n;
  314. complex *zNew;
  315. {
  316.     complex quotient;
  317.     Divide(poly[n], poly[n-1], "ient);
  318.     Subtract(zOld, quotient, zNew);
  319. }
  320.  
  321.     
  322. Guess(z)
  323. /* A random complex number, 50% chance of being negative. */
  324. complex *z;
  325. {
  326. #ifdef UNIX
  327. #define    ran    drand48
  328. #endif
  329.     double ran();
  330.  
  331.     z->n[REAL] = ran();
  332.     z->n[IMAG] = ran();
  333.     z->n[REAL] = (ran() < 0.50) ? z->n[REAL] : -(z->n[REAL]);
  334.     z->n[IMAG] = (ran() < 0.50) ? z->n[IMAG] : -(z->n[IMAG]);
  335. }
  336.  
  337.     
  338. BOOL ErrorTooBig(y, z, eps)
  339. /* Compute the Euclidean distance between y and z in the complex plane.
  340.  * This is just our good old friend, the "distance formula".  (Add the
  341.  * squares of y and z, and take the square root.)  Is the distance less
  342.  * than epsilon?  If so, the error is allowable; else, it's too big. */
  343. complex y, z;
  344. double eps;
  345. {
  346.     complex diff;
  347.     double sqrt();
  348.  
  349.     Subtract(y, z, &diff);
  350.     return(    sqrt(SQR(diff.n[REAL]) + SQR(diff.n[IMAG])) < eps
  351.         ? FALSE
  352.         : TRUE);
  353. }
  354. \Rogue\Monster\
  355. else
  356.   echo "will not over write complex.c"
  357. fi
  358. if [ `wc -c complex.c | awk '{printf $1}'` -ne 2539 ]
  359. then
  360. echo `wc -c complex.c | awk '{print "Got " $1 ", Expected " 2539}'`
  361. fi
  362. if `test ! -s decl.h`
  363. then
  364. echo "writing decl.h"
  365. cat > decl.h << '\Rogue\Monster\'
  366. /* decl.h:    Declarations for the "Newton" program.
  367.  *
  368.  * Written by Daniel Barrett.  100% PUBLIC DOMAIN. */
  369.  
  370. #include <stdio.h>
  371. #include <math.h>
  372.  
  373. #ifdef UNIX
  374. #define    BOOL    short
  375. #define    TRUE    1
  376. #define    FALSE    0
  377. #else
  378. #include <exec/types.h>
  379. #endif
  380.     
  381. #define MAX_DEGREE    20        /* Maximum degree of equation. */
  382. #define    MAX_ITERATIONS    10000
  383. #define    EPSILON_DEFAULT    0.0000001
  384.  
  385. #define    EQUAL        !strcmp
  386. #define    SQR(x)        ((x)*(x))
  387.     
  388. struct ComplexNumber {            /* One complex number.     */
  389.     double n[2];            /* Real & imaginary parts. */
  390. };
  391. typedef struct ComplexNumber complex;
  392.  
  393. #define    REAL        0        /* 1st element of ComplexNumber. */
  394. #define    IMAG        1        /* 2nd element of ComplexNumber. */
  395.  
  396.  
  397. /* If C is of type "complex", then:
  398.  *
  399.  *    C.n[REAL]    is the real part.
  400.  *    C.n[IMAG]    is the imaginary part.
  401.  *
  402. */
  403. \Rogue\Monster\
  404. else
  405.   echo "will not over write decl.h"
  406. fi
  407. if [ `wc -c decl.h | awk '{printf $1}'` -ne 782 ]
  408. then
  409. echo `wc -c decl.h | awk '{print "Got " $1 ", Expected " 782}'`
  410. fi
  411. if `test ! -s ds.c`
  412. then
  413. echo "writing ds.c"
  414. cat > ds.c << '\Rogue\Monster\'
  415. /* ds.c:    Data structure routines.
  416.  *
  417.  * Written by Daniel Barrett.  100% PUBLIC DOMAIN. */
  418.  
  419. #include "decl.h"
  420.  
  421.  
  422. InitPoly(poly, n)
  423. /* Initialize the polynomial to all zeroes. */
  424. complex poly[];
  425. int n;
  426. {
  427.     int i;
  428.     for (i=0; i<n; i++)
  429.         AssignComplex(&poly[i], 0.0, 0.0);
  430. }
  431.  
  432.  
  433. AssignComplex(comp, realPart, imagPart)
  434. /* Assign the real & imaginary parts to a complex number. */
  435. complex *comp;
  436. double realPart, imagPart;
  437. {
  438.     comp->n[REAL] = realPart;
  439.     comp->n[IMAG] = imagPart;
  440. }
  441.  
  442.  
  443. CopyPoly(poly1, poly2, degree)
  444. /* Copy poly1 into poly2. */
  445. complex poly1[], poly2[];
  446. int degree;
  447. {
  448.     int i;
  449.     for (i=0; i<=degree; i++)
  450.         CopyComplex(poly1[i], &poly2[i]);
  451.  
  452. }
  453.  
  454.  
  455. CopyComplex(c1, c2)
  456. /* Copy complex number c1 into c2. */
  457. complex c1, *c2;
  458. {
  459.     AssignComplex(c2, c1.n[REAL], c1.n[IMAG]);
  460. }
  461.  
  462. \Rogue\Monster\
  463. else
  464.   echo "will not over write ds.c"
  465. fi
  466. if [ `wc -c ds.c | awk '{printf $1}'` -ne 775 ]
  467. then
  468. echo `wc -c ds.c | awk '{print "Got " $1 ", Expected " 775}'`
  469. fi
  470. if `test ! -s io.c`
  471. then
  472. echo "writing io.c"
  473. cat > io.c << '\Rogue\Monster\'
  474. /* io.c:    Input/Output functions.
  475.  *
  476.  * Written by Daniel Barrett.  100% PUBLIC DOMAIN. */
  477.  
  478. #include "decl.h"
  479.  
  480.  
  481. PrintComplex(comp)
  482. /* Print a complex number. */
  483. complex comp;
  484. {
  485.     double fabs();
  486.  
  487.     if (comp.n[REAL] >= 0)
  488.         printf(" ");
  489.     printf("%f%c%fi",    comp.n[REAL], 
  490.                 (comp.n[IMAG] < 0.0) ? '-' : '+',
  491.                 fabs(comp.n[IMAG]));
  492. }
  493.  
  494.     
  495. PrintPoly(poly, n)
  496. /* Print the polynomial in a nice format. */
  497. complex poly[];
  498. int n;
  499. {
  500.     int i;
  501.  
  502.     printf("\nYour polynomial is:\n");
  503.     printf("--------------------\n");
  504.     for (i=0; i<n-1; i++) {
  505.         PrintComplex(poly[i]);
  506.         printf("    * X^%d    +\n", n-i);
  507.     }
  508.     PrintComplex(poly[n-1]);
  509.     printf("    * X    +\n");
  510.     PrintComplex(poly[n]);
  511.     printf("\n\n");
  512. }
  513.  
  514.     
  515. ReadPoly(poly, degree, epsilon)
  516. /* Read all data from the user:  the degree of the polynomial, the
  517.  * accuracy desired, and the polynomial coefficients. */
  518. complex poly[];
  519. int *degree;
  520. double *epsilon;
  521. {
  522.     ReadDegree(degree);
  523.     ReadEpsilon(epsilon);
  524.     TellUserWhatToDo();
  525.     ReadCoefficients(poly, *degree);
  526. }
  527.  
  528.     
  529. ReadDegree(degree)
  530. /* Prompt user for the degree of the polynomial.  Read it. */
  531. int *degree;
  532. {
  533.     char buf[BUFSIZ];
  534.  
  535.     *degree = 0;
  536.     do {
  537.         printf("Degree of polynomial? (1..%d): ", MAX_DEGREE);
  538.         gets(buf);
  539.         *degree = atoi(buf);
  540.     }while (*degree <= 0 || *degree > MAX_DEGREE);
  541. }
  542.  
  543.     
  544. ReadEpsilon(epsilon)
  545. /* Prompt user for epsilon, the desired accuracy.  Read it.  If no
  546.  * data, set epsilon to the default value. */
  547. double *epsilon;
  548. {
  549.     char buf[BUFSIZ];
  550.     double atof();
  551.  
  552.     do {
  553.         *epsilon = EPSILON_DEFAULT;
  554.         printf("Accuracy desired? [<RETURN> = %1.7f]: ", *epsilon);
  555.         gets(buf);
  556.          if (strlen(buf))
  557.             *epsilon = atof(buf);
  558.     }while (*epsilon <= 0.0);
  559. }
  560.  
  561.     
  562. ReadCoefficients(poly, degree)
  563. /* Read in the coefficients of the polynomial "poly".  NOTE that
  564.  *  the array "poly" is backwards... poly[0] is the coefficient
  565.  *  of X^degree, and so on. 
  566.  * I use "strtok()", my version of the UNIX function, to get the input. */
  567. complex poly[];
  568. int degree;
  569. {
  570.     char buf[BUFSIZ], *tokeReal=NULL, *tokeImag=NULL, *strtok();
  571.     static char separators[] = " \t\n\r";
  572.     register int i;
  573.     double atof();
  574.  
  575.     for (i=0; i <= degree; i++) {
  576. Again:        printf("X^%d coefficient: ", degree-i);
  577.         gets(buf);
  578.         tokeReal = strtok(buf, separators);
  579.         if (!tokeReal)
  580.             goto Again;
  581.         else if (NeedsHelp(tokeReal)) {
  582.             Help();
  583.             goto Again;
  584.         }
  585.         tokeImag = strtok((char *)NULL, separators);
  586.         AssignComplex(&poly[i],    atof(tokeReal),
  587.                     tokeImag ? atof(tokeImag) : 0.0);
  588.     }
  589. }
  590.  
  591.  
  592. PrintRoot(root, i)
  593. /* Print the value of root "i" of the polynomial. */
  594. complex root;
  595. int i;
  596. {
  597.     printf("Root %d:  ", i);
  598.     if (i<10)
  599.         printf(" ");
  600.     PrintComplex(root);
  601.     printf("\n");
  602. }
  603.  
  604.     
  605. TellUserWhatToDo()
  606. {
  607.     printf("Enter the coefficients of your polynomial.  ");
  608.     printf("Type \"H\" for help.\n");
  609. }
  610.  
  611.     
  612. static char *helpMessage[] = {
  613. "",
  614. "Enter the coefficients of your polynomial, one coefficient per line.",
  615. "Note that \"X^n\" means \"X to the nth power\".",
  616. "If the coefficient is REAL, simply enter it and press <RETURN>.",
  617. "If the coefficient is COMPLEX, enter its real and imaginary terms,",
  618. " separated by a space.  (Do not use the letter \"i\".)",
  619. "",
  620. NULL
  621. };
  622.  
  623.     
  624. Help()
  625. /* Print the above help message. */
  626. {
  627.     char **s = helpMessage;
  628.     do
  629.         puts(*s);
  630.     while (*(s++));
  631. }
  632.  
  633.  
  634. NeedsHelp(s)
  635. /* Does the user need help? */
  636. char *s;
  637. {
  638.     return(EQUAL(s, "H") | EQUAL(s, "h") | EQUAL(s, "?"));
  639. }
  640.  
  641.     
  642. Version()
  643. {
  644.     printf("[33mNewton V1.0[0m");
  645.     printf(" by Daniel Barrett.  100%% PUBLIC DOMAIN.\n");
  646.     printf("Estimate the roots of a polynomial numerically.\n");
  647.     printf("Type ^C <RETURN> to abort this program.\n\n");
  648. }
  649. \Rogue\Monster\
  650. else
  651.   echo "will not over write io.c"
  652. fi
  653. if [ `wc -c io.c | awk '{printf $1}'` -ne 3554 ]
  654. then
  655. echo `wc -c io.c | awk '{print "Got " $1 ", Expected " 3554}'`
  656. fi
  657. if `test ! -s main.c`
  658. then
  659. echo "writing main.c"
  660. cat > main.c << '\Rogue\Monster\'
  661. /* main.c:    Main "Newton" program.
  662.  *
  663.  * Written by Daniel Barrett.  100% PUBLIC DOMAIN. */
  664.  
  665. #include "decl.h"
  666.     
  667. main(argc, argv)
  668. int argc; char *argv[];
  669. {
  670.     complex polynomial[MAX_DEGREE+1];    /* Our polynomial.   */
  671.     int degree;                /* Its degree.       */
  672.     double epsilon;                /* Desired accuracy. */
  673.  
  674. #ifdef UNIX
  675.     srand48((long)123);
  676. #endif
  677.  
  678.  
  679.     epsilon    = 0.0;
  680.     degree    = 0;
  681.  
  682.     Version();
  683.     InitPoly(polynomial, MAX_DEGREE+1);
  684.     ReadPoly(polynomial, °ree, &epsilon);
  685.     PrintPoly(polynomial, degree);
  686.     NewtonMethod(polynomial, degree, epsilon);
  687. }
  688.  
  689.  
  690. NewtonMethod(poly, degree, epsilon)
  691. /* Find the roots using Newton's Method. */
  692. complex poly[];
  693. int degree;
  694. double epsilon;
  695. {
  696.     int numRoots=0, iterations;
  697.     complex zero, oldGuess, newGuess, polyCopy[MAX_DEGREE+1];
  698.     BOOL ErrorTooBig();
  699.  
  700.     AssignComplex(&zero, 0.0, 0.0);
  701.     while (numRoots < degree-1) {        /* For each root...   */
  702.         CopyPoly(poly, polyCopy, degree);
  703.         Guess(&oldGuess);        /* Guess its value... */
  704.         CopyComplex(oldGuess, &newGuess);
  705.         iterations = MAX_ITERATIONS;
  706.         do {                /* Refine the guess...*/
  707.             CopyPoly(poly, polyCopy, degree);
  708.             CopyComplex(newGuess, &oldGuess);
  709.             SynthDivide(polyCopy, oldGuess, degree-numRoots);
  710.             SynthDivide(polyCopy, oldGuess, (degree-1)-numRoots);
  711.             Iterate(polyCopy, oldGuess, degree-numRoots, &newGuess);
  712.         }while (ErrorTooBig(oldGuess, newGuess, epsilon)
  713.             &&  --iterations);        /* ...until convergence. */
  714.         if (!iterations) {
  715.             printf("Root didn't converge after %d iterations.\n",
  716.                 MAX_ITERATIONS);
  717.             printf("Exiting!\n");
  718.             exit(20);
  719.         }
  720.         SynthDivide(poly, newGuess, degree-numRoots);
  721.         PrintRoot(newGuess, ++numRoots);
  722.     }
  723.  
  724.     Divide(poly[1], poly[0], &poly[1]);    /* Get the last root. */
  725.     Subtract(zero, poly[1], &newGuess);
  726.     PrintRoot(newGuess, degree);
  727. }
  728. \Rogue\Monster\
  729. else
  730.   echo "will not over write main.c"
  731. fi
  732. if [ `wc -c main.c | awk '{printf $1}'` -ne 1760 ]
  733. then
  734. echo `wc -c main.c | awk '{print "Got " $1 ", Expected " 1760}'`
  735. fi
  736. if `test ! -s strtok.c`
  737. then
  738. echo "writing strtok.c"
  739. cat > strtok.c << '\Rogue\Monster\'
  740. /* strtok.c:    3 very useful string-handling functions.
  741.  *        You don't need this file if you are using Lattice C.
  742.  *        Manx C does not provide them.
  743.  *
  744.  * Written by Daniel Barrett.  100% PUBLIC DOMAIN. */
  745.  
  746.  
  747. #include "decl.h"
  748. #define STRING_END    '\0'
  749.  
  750. /* Return the next string "token" in "buf".  Tokens are substrings
  751.  * separated by any characters in "separators". */
  752.  
  753. char *strtok(buf, separators)
  754. char *buf, *separators;
  755. {
  756.     register char *token, *end;    /* Start and end of token. */
  757.     extern char *strpbrk();
  758.     static char *fromLastTime;
  759.  
  760.     if (token = buf ? buf : fromLastTime) {
  761.         token += strspn(token, separators);    /* Find token! */
  762.         if (*token == STRING_END)
  763.             return(NULL);
  764.         fromLastTime = ((end = strpbrk(token,separators))
  765.                 ? &end[1]
  766.                 : NULL);
  767.         *end = STRING_END;            /* Cut it short! */
  768.     }
  769.     return(token);
  770. }
  771.  
  772.     
  773. /* Return a pointer to the first occurance in "str" of any character
  774.  * in "charset". */
  775.  
  776. char *strpbrk(str, charset)
  777. register char *str, *charset;
  778. {
  779.     extern char *index();
  780.  
  781.     while (*str && (!index(charset, *str)))
  782.         str++;
  783.     return(*str ? str : NULL);
  784. }
  785.  
  786.     
  787. /* Return the number of characters from "charset" that are at the BEGINNING
  788.  * of string "str".
  789. */
  790.  
  791. int strspn(str, charset)
  792. register char *str, *charset;
  793. {
  794.     register char *s;
  795.     s = str;
  796.     while (index(charset, *s))
  797.         s++;
  798.     return(s - str);
  799. }
  800.  
  801. \Rogue\Monster\
  802. else
  803.   echo "will not over write strtok.c"
  804. fi
  805. if [ `wc -c strtok.c | awk '{printf $1}'` -ne 1322 ]
  806. then
  807. echo `wc -c strtok.c | awk '{print "Got " $1 ", Expected " 1322}'`
  808. fi
  809. echo "Finished archive 1 of 1"
  810. # if you want to concatenate archives, remove anything after this line
  811. exit
  812. -- 
  813. Bob Page, U of Lowell CS Dept.  page@swan.ulowell.edu  ulowell!page
  814. Have five nice days.
  815.